
**********
* Readme *
**********

* This script 
* [1] estimates a quantile regression model (5th, 10th 15ht percentiles) on pof data and
* [2] fits the expected reservation wage conditional on individual's covariates.


* Root folder (PATH TO BE DEFINED BY THE USER)
**********************************************
clear all
global analysis "C:/***/replication_package"


* Timestamped log
*****************
global today = strofreal(date(c(current_date), "DMY"), "%tdYYNNDD")
log using "${analysis}/code/logs/4_2_estimate_reservation_wages_${today}.smcl", replace


**************
* Import POF *
**************

use "${analysis}/data/2_2_pof_clean.dta", clear


* Covariates
************

* Key idea: what's the lowest wages usually accepted by people similar to me?

* Define set of covariates
unab covariates_all : rg_* ae_* at_* fp_* n_kids n_youngs n_adults n_seniors region_* // all
unab ref_groups : rg_1_nonwhite_female ae_1_1 fp_h_wp_nk region_SP_c // reference group: non white female, no education, head w partner no kids, Sao Paulo
global covariates : list covariates_all - ref_groups // final list of regressors


*******************************************
* Baseline estimation & alternative specs *
*******************************************

* Preferred definition (domestic workers as OAWs)
replace ln_winc_employee = .
replace ln_winc_employee = ln_winc if work_state_name == "Employee"

qreg ln_winc_employee $covariates [pweight = pweight], quantile(.05) wlsiter(50) iterate(600) vce(robust, chamberlain)
estimates save "${analysis}/results/estimations/4_2_mod_reserv_wage_05.ster", replace

qreg ln_winc_employee $covariates [pweight = pweight], quantile(.10) wlsiter(50) iterate(600) vce(robust, chamberlain)
estimates save "${analysis}/results/estimations/4_2_mod_reserv_wage_10.ster", replace

qreg ln_winc_employee $covariates [pweight = pweight], quantile(.15) wlsiter(50) iterate(600) vce(robust, chamberlain)
estimates save "${analysis}/results/estimations/4_2_mod_reserv_wage_15.ster", replace

* IBGE default definition (domestic workers as employees)
replace ln_winc_employee = .
replace ln_winc_employee = ln_winc if work_state_ibge_name == "Employee"

qreg ln_winc_employee $covariates [pweight = pweight], quantile(.10) wlsiter(100) iterate (500) vce(robust, chamberlain)
estimates save "${analysis}/results/estimations/4_2_appx_dom_mod_reserv_wage_10.ster", replace


***********************************************
* Fitted values: baseline & alternative specs *
***********************************************

use "${analysis}/data/2_2_pof_clean.dta", clear
svyset psu_id [pweight = pweight], strata(strata_id) singleunit(centered) vce(linearized) 

estimates use "${analysis}/results/estimations/4_2_mod_reserv_wage_05.ster"
predict est_ln_rw_05_0, xb
predictnl est_rw_05_0 = exp(xb())

estimates use "${analysis}/results/estimations/4_2_mod_reserv_wage_10.ster"
predict est_ln_rw_10_0, xb
predictnl est_rw_10_0 = exp(xb())

estimates use "${analysis}/results/estimations/4_2_mod_reserv_wage_15.ster"
predict est_ln_rw_15_0, xb
predictnl est_rw_15_0 = exp(xb())

estimates use "${analysis}/results/estimations/4_2_appx_dom_mod_reserv_wage_10.ster"
predict appx_dom_est_ln_rw_10_0, xb
predictnl appx_dom_est_rw_10_0 = exp(xb())

label var          est_rw_05_0 "Estimated reserv wage (q = 5)"
label var          est_rw_10_0 "Estimated reserv wage (baseline)"
label var          est_rw_15_0 "Estimated reserv wage (q = 15)"
label var appx_dom_est_rw_10_0 "Estimated reserv wage (dom worker as employee)"

label var          est_ln_rw_05_0 "Estimated log reserv wage (q = 5)"
label var          est_ln_rw_10_0 "Estimated log reserv wage (baseline)"
label var          est_ln_rw_15_0 "Estimated log reserv wage (q = 15)"
label var appx_dom_est_ln_rw_10_0 "Estimated log reserv wage (dom worker as employee)"

* Summary os estimations
svy, subpop(if work_state_name == "Employee"):                  mean est_rw_05_0 est_rw_10_0 est_rw_15_0 appx_dom_est_rw_10_0
svy, subpop(if work_state_name == "Own-account worker"):        mean est_rw_05_0 est_rw_10_0 est_rw_15_0 appx_dom_est_rw_10_0
svy, subpop(if work_state_ibge_name == "Employee"):             mean est_rw_05_0 est_rw_10_0 est_rw_15_0 appx_dom_est_rw_10_0
svy, subpop(if work_state_ibge_name == "Own-account worker"):   mean est_rw_05_0 est_rw_10_0 est_rw_15_0 appx_dom_est_rw_10_0


* Store the fitted values
*************************
keep  ind_id est_ln_rw_* est_rw_* appx_dom_est_ln_rw_10_0 appx_dom_est_rw_10_0
order ind_id est_ln_rw_* est_rw_* appx_dom_est_ln_rw_10_0 appx_dom_est_rw_10_0
describe
save "${analysis}/data/4_2_est_reservation_wages.dta", replace


****************************
* Estimation output tables *
****************************

use "${analysis}/data/2_2_pof_clean.dta", clear
svyset psu_id [pweight = pweight], strata(strata_id) singleunit(centered) vce(linearized) 
keep if !missing(ln_winc_oaw)

estimates use "${analysis}/results/estimations/4_2_mod_reserv_wage_05.ster"
estimates store qreg_05

estimates use "${analysis}/results/estimations/4_2_mod_reserv_wage_10.ster"
estimates store qreg_10

estimates use "${analysis}/results/estimations/4_2_mod_reserv_wage_15.ster"
estimates store qreg_15

global region region_*

* Print
esttab qreg_05 qreg_10 qreg_15, varwidth(50) ///
  cells("b(fmt(3)) se(par fmt(3))") ///
  label mlabels(, titles) nonumbers eqlabels(none) collabel("coef." "s.e.") ///
  drop($region _cons) ///
  refcat( ///
    rg_2_white_female "Ethnicity and gender (ref: Nonwhite female)" ///
    ae_1_2 "Age and education (ref: 14-24, less than prim. school)" ///
    ae_2_1 "[BLANK LINE]" ///
    ae_3_1 "[BLANK LINE]" ///
    ae_4_1 "[BLANK LINE]" ///
    ae_5_1 "[BLANK LINE]" ///
    at_school "Current schooling status (ref: Not currently studying)" ///
    fp_h_wp_wk "Household position (ref: Head, with partner, no kids)" ///
    n_kids "Number of household members by age", nolabel)


* TABLE 6
*********

esttab qreg_05 qreg_10 qreg_15 using "${analysis}/results/tables/table6.tex", style(tex) fragment replace ///
  cells("b(fmt(3)) se(par fmt(3))") ///
  label mgroups(none) mlabels(none) nonumbers eqlabels(none) collabel(none) noobs ///
  drop($region _cons) ///
  refcat( ///
    rg_2_white_female   "[nospace]\textit{Ethnicity and gender} [emptycolumns] \\ \hspace{2ex}Nonwhite female" ///
    ae_1_2              "[nospace]\textit{Age and education} [emptycolumns] \\ \hspace{2ex}14-24, less than prim. school" ///
    at_school           "[nospace]\textit{Current schooling status} [emptycolumns] \\ \hspace{2ex}Not currently studying" ///
    fp_h_wp_wk          "[nospace]\textit{Household position} [emptycolumns] \\ \hspace{2ex}Head, with partner, no kids" ///
    n_kids              "[nospace]\textit{Number of household members by age}", nolabel) ///
  varlabels(, /// 
    prefix(\hspace{2ex}) ///
    elist(rg_4_white_male \tabularnewline ae_1_4 \tabularnewline ae_2_4 \tabularnewline ae_3_4 \tabularnewline ae_4_4 \tabularnewline ae_5_4 \tabularnewline at_college \tabularnewline fp_oa \tabularnewline)) ///
  substitute("\hline" "" ///
             "&" " & " "  " " " "  " " " "  " " " "  " " " ///
             " & & " " & $\cdot$ & " ///
             " & & " " & $\cdot$ & " ///
             " & & " " & $\cdot$ & " ///
             " & \\" " & $\cdot$ \\" ///
             "\hspace{2ex}[nospace]" "" ///
             "[emptycolumns]" "& & & & & &" ///
             "\textit{Number of household members by age} & $\cdot$ & $\cdot$ & $\cdot$ & $\cdot$ $\cdot$ & $\cdot$ & \\" "\textit{Number of household members by age} & & & & & & \\")

             
             
* End of script
***************
cap log close